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Recently, it has been suggested that the Many-Body Localized phase can be characterized by 
local integrals of motion. Here we introduce a Hilbert space preserving renormalization scheme that 
iteratively finds such integrals of motion exactly. Our method is based on the consecutive action of a 
similarity transformation using displacement operators. We show, as a proof of principle, localization 
and the delocalization transition in interacting fermion chains with random onsite potentials. Our 
scheme of consecutive displacement transformations can be used to study Many Body Localization 
in any dimension, as well as disorder-free Hamiltonians. 
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Since the revival of interest in localization due to 
disorder [IHl] it has been suggested that the so-called 
Many-Body Localized (MBL) phase can be characterized 
by an extensive set of local integrals of motion (LIOM) or 
Lbits, r/, that commute with each other and the Hamil¬ 
tonian EiiMn]. Consequently, the Hamiltonian can be 
written in terms of these LIOMs as 

H = + ( 1 ) 

i ij 

Many properties of the MBL phase, such as its logarith¬ 
mic entanglement spread or its insulating behavior, can 
be derived based on this assumption m- 

The question is, however, what those LIOMs are and 
how to compute them. Since any sum and product of in¬ 
tegrals of motion is itself an integral of motion, the choice 
of LIOMs is highly arbitrary. Pure mathematically, all 
projectors onto the (localized) eigenstates are integrals of 
motion, and out of those one could in principle construct 
the local integrals of motion. In fact, it is easy to show 
that all Hamiltonians can be brought into the form dic¬ 
tated by Eqn. 0 m- As for the MBL phase, Chandran 
et al.[3] use the long-time evolved average of an initially 
local operator as their LIOMs, whereas Ros et al.[S] and 
Imbrie[7] use perturbative methods to construct local in¬ 
tegrals of motion. 

In this Letter, we construct iteratively a transforma¬ 
tion that turns any fermionic Hamiltonian into the clas¬ 
sical form of Eqn. 0. This is done by consecutively ap¬ 
plying a similarity transformation using a displacement 
operator expA(Xl — X). The elegant properties of this 
transformation allow for a systematic elimination of off- 
diagonal interaction terms, order by order in the num¬ 
ber of fermionic operators involved. Our renormalization 
scheme can be used to study Hamiltonians in any dimen¬ 
sion, with or without disorder. 

Whether a random interacting system is localized or 
not, depends on how much the integrals of motion rf are 
spread out. As a proof of principle, we apply our method 
to diagonalize random interacting chains. In both the lo¬ 


calized and the delocalized regime we find agreement with 
Exact Diagonalization. Throughout the phase diagram 
we find the effective interactions between the integrals of 
motion, and we can infer the exponential localization of 
the integrals of motion in the localized regime. 

Definitions - Here we will consider interacting 
fermions with random onsite potentials [2113] 

H — ^ ^ + 2 ^ ^ Ed/S7(5cjjC^C-yC5, (2) 

a a/S-yS 

as it extends the original concept of Anderson 
localization [T] to interacting systems. For this model we 
will now characterize all the possible terms in the Hamil¬ 
tonian. 

We define a classical term in the Hamiltonian as 
a product of fermionic density operators of the form 
ni.,^ni^ .. .Ui^ where all i's are different. A quantum term 
is the product of fermionic operators that cannot be writ¬ 
ten as a classical term. It thus contains, next to a possible 
set of density operators, separate creation and annihila¬ 
tion operators, 

X = ...Cj^. ( 3 ) 

Again, all the i’s and all the j’s are different from each 
other. The order 0{X) of a term is defined as the to¬ 
tal number of creation and annihilation operators it con¬ 
tains. This equals 0{X) = 2d + q where d is the num¬ 
ber of density operators and q the number of separate 
creation and annihilation operators. We require the or¬ 
der to be even, and the interaction to contain the same 
number of creation and annihilation operators, so that all 
interactions preserve the total fermion number |21|. We 
call the notation of Eqn. ([^ normal ordered: all density 
interactions are grouped together, and the remaining is 
a product of alternating creation and annihilation oper¬ 
ators which all act on different sites. 

The order of a sum of terms is defined as the minimum 


2 


of the orders of the individual terms, 

O =min(0(X,)). (4) 

When multiplying two terms X and Y, the product 
XY can contain interaction terms of order lower than 
0{X) + 0{Y). This happens when X contains the anni¬ 
hilation operator on site a and Y contains the creation 
operator on the same site, we call this an overlap. This 
gives Cac\ = 1 — ria, which generates a term of an order 
two lower. Since this can happen for any pair of creation 
and annihilation operators in X and Y, there exists a case 
of maximal overlap with overlap on 4 min(C>(X), 0(y)) 
sites. Because each such overlap generates a term with 
an order two lower, the product XY has order 

0{X) + 0{Y) > 0{XY) > max {0{X),0{Y)). (5) 

There are three important properties of quantum 
terms: 1) The product of a quantum term with it¬ 
self is zero, X^X'^ = XX = 0. 2) The product of a 
quantum term with its Hermitian conjugate, so X^X 
and XX^, is classical. Observe that because X and 
X^ have maximal overlap, the order remains the same: 
OiXX-i) = 0{X^X) = OiX). 3) The cubic power is 
trivial, that is X^XX^ = X^ and XX^ X = X. 

Displacement Transformation - Based on the classifi¬ 
cation of terms we just introduced, we can define a dis¬ 
placement operator associated with a quantum term X, 

Dx{X)=exp{X{X^ -X)). (6) 

Because of the aforementioned properties of quantum 
terms, the displacement operator can be written out ex¬ 
plicitly as 

Vx{X) = 1 + sinXiX'^ - X) + {cos X - 1){X^ X + XX^). 

(7) 

Note that the Hermitian conjugate of the displacement 
operator is dI^{X) = Dx{—X). The above arguments 
can easily be extended to spin-4 Hamiltonians, where 
classical terms are given by products of 5'^-operators, and 
quantum terms are total spin-conserving products of S~^ , 
S~ and operators. 

The displacement transformation is given by a sim¬ 
ilarity transformation using the displacement operator. 
That is, it transforms any term Y as 

Y ^Y = vl,iX)YVx{X). (8) 

This transformation is similar to a Clifford group rota¬ 
tion, which has transformation operator D = with 
= 1, whereas here we have the weaker condition 
= -A. 

For now we use the notation with the tilde to denote 
the transformed term, we will drop the tilde later as 


we will perform many consecutive transformations. Un¬ 
der the transformation, there are ’new’ terms generated^ 
namely Y — Y. Using the explicit formulation of the dis¬ 
placement operator Eqn. ([^ , we see that the ’new’ terms 
are of the form XY, YX, etc. Carefully counting all the 
combinations, we see that the order of the new terms is 
at least the maximum of the orders of X and Y, 

C>(y-y) >max(C>(X),C>(y)). (9) 

As will be shown later, this lower bound on the order 
of new terms implies the closedness of our systematic 
transformation procedure. 

Without constraining the specific shape of X that we 
use for the displacement transformation, we can prove 
that the only way to generate new terms proportional 
to X^ + X is through terms that have maximal over¬ 
lap with xm- For example, consider an order 4 term 
At = d C 2 C 3 C 4 with interaction given by V. Let us write 
down the relevant part of the Hamiltonian as follows 

4 1 

-t- Visuius -t- V24n2n4 -t- ( 1 ^) 

i=l 


The displacement transformation with Dx{X) leaves the 
quadratic part untouched, and the prefactor multiplying 
(At -I- A) becomes 

-V cos 2A -I- - (^1 -I- ^3 -f V 13 — ^2 ~ ^4 ~ 1 ^ 24 ) sin 2A 

( 11 ) 

so that with A given by 


tan2A = ----—-----. 

Cl + ^3 + ^13 — C2 — Cl — V24 


( 12 ) 


the transformed Hamiltonian does no longer have the in¬ 
teraction term At X. Similar expressions can be found 
for transformations involving A of higher order. 

The right-hand side of Eqn. (12) equals the ’small’ pa¬ 
rameter that is used in perturbative studies of MBL[^|31 
[7] . Such perturbation theories often run into the problem 
of resonances, where the denominator of the ’small’ pa¬ 
rameter goes to zero, which means perturbation theory 
cannot be applied. However, the displacement transfor¬ 
mation we present here is well-behaved at a resonance, 
since then A = 7 r /4 and the interaction can be still trans¬ 
formed away. Note that each displacement transforma¬ 
tion can be viewed as a discrete version of a Wegner 
transformation [52] . 

Consecutive displacement transformations - Now any 
fermionic Hamiltonian that respects the total fermion 
number conservation can be written as 


= ^ Y.^uj{xl^ + X^,) (13) 

i n—4,6... j 


where Ci are the onsite energies, n expresses the order of 
the term A„j, j is just an index and U„j are the coupling 
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constants. If all terms X are classical, we have reached 
our goal: we have a classical Hamiltonian with an infinite 
set of conserved quantities. 

Any quantum term + A") can be removed from 
the Hamiltonian by performing a displacement transfor¬ 
mation associated with X, using the value of A given by 


Eqn. (121. After done so, we can choose another quan¬ 


tum term and transform that one away - and continue 
this path of consecutive transformations. 

New terms that are generated are multiplied by either 
sin(A), cos A or products of those. Therefore, generically, 
new terms have smaller couplings constants, making the 
process of consecutive transformations alike a renormal¬ 
ization scheme. This generation of new terms takes into 
account co-tunneling and hence long distance resonances. 
Note that our scheme preserves the Hilbert space, since 
each displacement transformation does not decimate sites 
or bonds nor any other degree of freedom. 

In certain cases, however, transforming a term X away 
can generate terms with an even larger coupling constant. 
This does not pose a problem: whenever we transform a 
term X away and it is later regenerated, upon regener¬ 
ation it will have a smaller coupling constant than be¬ 
fore. Additionally, throughout consecutive transforma¬ 
tions the distribution of coupling constants will change 
such that creating larger coupling constants will become 
more and more unlikely. Therefore the magnitude of the 
strongest coupling constant decreases exponentially with 
the number of applied transformations |21) . 

We thus remove, term by term, all quantum terms of 
order 4 in the Hamiltonian. The price we pay is the gener¬ 
ation of new terms (both classical and quantum) of order 
6 and higher, and new classical terms of order 4. As a 
result, we obtain a complicated Hamiltonian that is clas¬ 
sical in its quadratic and quartic terms. Subsequently, 
we can do the same tricks for the next order quantum 
terms, making the Hamiltonian at that level classical as 
well, and so forth. The procedure cut off at n-th order re¬ 
produces the exact spectrum for states with n/2 particles 
or less. The computational complexity of diagonalizing 
a Hamiltonian up to n-th order is ©(A^^"/^), which we 
confirmed numerically [H]. 

In practice one needs to cut off the procedure at a 
certain order n. This approximation has a clear physi¬ 
cal interpretation: we are expressing many-body states 
systematically in terms of n-particle states. 

Numerical implementation - As a proof of principle, 
we implemented our method numerically. We consider 
an open chain of L sites with spinless fermions, with a 
random onsite energy Ci on each site chosen uniformly 
between —Wjl and IH/2, hopping t = 1 and a nearest 
neighbor repulsion with V = 1, 


N 


N-1 




iV-l 




( 14 ) 



FIG. 1: (Color online) The average of the absolute value of 
the coefficients (| |) present in the classical Hamiltonian. 

Here we show results deep in the many-body localized phase 
with IT = 8, for a system size L = 24 averaged over 20 disor¬ 
der realizations. The different curves represent the two-body 
Vij, three-body Vijk and four-body Vijki coefficients, respec¬ 
tively. The strength of the coupling constants is almost inde¬ 
pendent of the order of the interaction. Inset: The spread of 
the local integrals of motion at disorder strength W — S and 
L = 20, compared to the coefficients of the Hamiltonian at the 
same distance. Note the difference between the localization 
length and the exponential decay of the interactions. 


The Hamiltonian is first diagonalized at the quadratic 
level, after that we continue with displacement transfor¬ 
mations at quartic order. 

At each iteration we pick the quantum term with the 
largest coupling constant at a given order, and transform 
it away. We neglect coupling constants smaller than nu¬ 
merical accuracy, set at e = 3 x 10“^. The Hamiltonian is 
thus diagonalized order by order, up to 8th order, yield¬ 
ing 

H = + + + (15) 

The results are averaged over 20 disorder realizations for 
L = 24 length chains, and 30 realizations at L = 20. 

Deep in the disordered phase, for IT = 8, the coefh- 
cients in the Hamiltonian fall off exponentially with the 
maximum distance of the sites involved. This is quan¬ 
tified by computing the average of the absolute value of 
the coefficients, {\Vi^,,,i„\) as shown in Fig. Note that 
the strength of the coupling constants appears to be in¬ 
dependent of the number of density terms involved. 

One advantage of our method is that we can also com¬ 
pute the classical Hamiltonian in the delocalized phase. 
In Fig. [^we show the distance-dependence of the Hamil¬ 
tonian coefficients as a function of disorder, ranging from 
IT = 9 to IT = 3. If the disorder is weaker than HA « 5 
the coefficients become independent of distance, which 
signals delocalization. 
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FIG. 2: (Color online) The average of the absolute value of 
the coefficients (| |) as a function of distance for various 

disorder strengths W = 3, ... ,9 for L = 20. In the localized 
phase one can see exponential decay, whereas in the delocal¬ 
ized phase at VF < Wc ~ 5 the coefficients are independent 
of distance. Inset: The relative error in the ground state en¬ 
ergy obtained using displacement transformations, compared 
to exact diagonalization results for L = 14. The different 
curves represent the cut-off at a given order n = 4,6, 8. In 
the localized phase even few-body interactions are sufficient 
to reproduce the ED results. In the delocalized regime one 
needs to go up to order 8, whereas around the MBL transition 
interactions of even higher order need to be included. 

To test the accuracy of our method, we compared the 
ground state energy obtained using our method with Ex¬ 
act Diagonalization (ED) results for L = 14, see the inset 
of Fig. With only terms up to 4th order we get a more 
than 95% accuracy in the localized phase. However, close 
to the transition and in the delocalized regime higher or¬ 
der terms are necessary to approach the ED results. For 
the delocalized regime the 8 th order terms seem to be suf¬ 
ficient, yet at the critical point we need even more-body 
interactions. This suggests our method becomes very ex¬ 
pensive close to the localization-delocalization transition. 

Having thus established that the practical implemen¬ 
tation of our scheme correctly reproduces ED results in 
the localized regime, we can directly probe the locality 
of the integrals of motion. The integrals of motion are 
given by transforming the initial density-operators 

rf = C/'l'n.t/ = n^+ <^)kifn c]ckc\c„, + . (16) 

jklm 

where U is the product of all the displacement transfor¬ 
mations. Note that at quadratic order, rf = rii, con¬ 
sistent with the fact that the single-particle spectrum is 
unchanged by the interactions. 

The average of the absolute value of the coefficients 
q; as a function of the distance away from the original 
site i is shown in the inset of Fig. The integrals of 
motion are indeed exponentially localized, as expected. 
The localization length for is, however, different from 


the length associated with the decay of the interactions 
coefficients Vij. 

These results serve as a proof of principle that the 
method of displacement transformations can be used to 
study interacting fermion models. In the Supplementary 
Information we also applied our method to a different 
modeipT]. In future work we will apply this method to 
extract new physical results for various models. 

Outlook - We introduced a sequence of displacement 
transformations that allows for the diagonalization of an 
interacting fermionic Hamiltonian. This method suggests 
we can bring any charge-conserving Hamiltonian into the 
classical form of Eqn. Q , not limited to the many-body 
localized phase[TT] or d = 1 . For example, even the com¬ 
pletely nonlocal Fermi liquids [13] can be analyzed using 
the classical model [T3j 

E = '^^knk + ^'^fkk'nknk' + ... (17) 

k kk' 

where the integrals of motion Uk are now local in momen¬ 
tum space and thus delocalized in real space. This begs 
the question how our method is related to the notion of 
integrability. We propose that in Eqn. ([^ quantum in- 
tegrability depends on the number of parameters 
that are nonzero, following the definition of Ref. [15]: if a 
subexponential or less number of parameters are nonzero 
and independent, the system is integrable. Note that 
this case of sparse uncorrelated coupling constants also 
implies Poissonian level statistics[T2]. Our method thus 
most effectively applies to systems close to integrability. 
It remains an open interesting question how ergodic sys¬ 
tems are represented in the r-basis of Eqn. 0. 

A possible fruitful future endeavor would be to recast 
the iterative transformations in the language of an ana¬ 
lytic renormalization scheme, much like the strong dis¬ 
order renormalization group theory [Tll 2 ], for which a re¬ 
lated version has been constructed for MBL systems [TU 
|5D]. A similar method using a Hilbert space preserv¬ 
ing RG scheme has been introduced by You et aip5]. 
The difference is that they keep interaction terms of all 
orders, and instead treat the off-diagonal resonance per- 
turbatively to second order. 

Even though the final classical Hamiltonian is of a re¬ 
markable simplicity, it does not imply easy solutions since 
in principle there could be long-range effective interac¬ 
tions between the classical bits r^. Furthermore, some 
complexity of the initial Hamiltonian is transferred to 
the transformation operator U, which is needed to trans¬ 
late any physical operator into the r-basis. Yet the fact 
that we can explicitly derive U using the renormalization 
scheme described in this Letter, introduces a novel quan¬ 
titative tool for the study of strongly interacting quantum 
matter. 
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MAXIMAL OVERLAP WITH THE TRANSFORMATION 

When a displacement transformation is performed with quantum term A, only terms with maximal overlap can 
generate new + X terms. Explicitly, we can enumerate all possible cases: 

1. A single density term rii where the site i corresponds to a creation operator in transforms as 

rii —>■ ni + i sin2A(A'f + A) — sin^ A(A'^A — AA'^). (18) 

The product of density terms nq ■ • • where all the sites correspond to creation operators in A^, 

generate the exact same new terms as the single density term. To prove this, observe that UiX^ = X^ and 
X^m = 0 . 

2. The same results hold, with —A instead of A, for a single density term Ui where the site i corresponds to a 
annihilation operator in A^ transforms as 

sin2A(Al' + A) + sin^ A(a1'A - AA^^). (19) 

This extends to products of density terms where the sites correspond to only annihilation operators in A^. 

3. The Hermitian interaction A^ + A transforms as 

(A^ + A) ^ cos2A(a 1' + A) -sin2A(Al'A - AAI'). (20) 

4. Under the displacement transformation with A, new terms of the form A^ + A can only be generated by the 
terms of the shape described in the previous three points. This follows from the fact that in order to generate 
terms like A^ + A, you need to transform a term Y that has order less than or equal to A, 0(Y) < 0{X), 
and which has maximal overlap with A. This means that all creation operators or all annihilation operators 
in Y should correspond to creation (annihilation) operators present in A or Xl Additionally, there cannot 
be operators in Y on sites that are not present in A, because those would be unaffected by the transition and 
all new terms would contain operators on these sites. Since Y also contains annihilation operators, they must 
either live on the same site as the creation operators (density terms), or they live on sites present in A. We 
have thus reduced the possible set of A’s to the three cases presented above. 


DIMINISHING OF LARGEST COUPLING CONSTANT 

In certain cases transforming a term A away can generate terms with an even larger coupling constant. We will 
now show that this does not cause a problem. Explicitly, imagine a Hamiltonian of the form 

= ^VxiX^ + A) + ^Vy{Y^ + A) + ^Uz(Zt + Z) + ... (21) 

such that Vx is the largest coupling constant, Vy, Vz < Vx- The terms Y and Z have a maximal overlap with A, 
and maximal overlap with each other. Upon transforming the term A away, the Hamiltonian becomes 

^ (Vy cos Ai + Uz sin Ai) (W^ + A) 

+ ^ {Vz cos Ai — Vy sin Ai) {Z^ + Z) + .. . (22) 

and it is clear that one (not both!) of the new coupling constants can be larger than the original Vx- If both coupling 
constants are smaller, we are contently moving closer to desired convergence. Instead, consider the unfortunate case, 
where VycosA + VysinA > Vx- By virtue of our system of consecutive transformations, the next step should be to 
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transform away the term Y. Doing so regenerates the original term X, however, this time with a smaller coupling 
constant, 

^ (Vz cos Ai — Vv sin Ai) sin A 2 (X'^ + X) 

+ ^ (VzcosAi — Vv sin Ai) cos A 2 (^^ + Z) + ... 

because (VzCOsAi — tVsinAi)sinA 2 < Vx by construction. This implies that the little detour caused by the larger 
coupling constant has come to an end, and the coupling constant in front of X has been reduced. 

Now one can track the magnitude of the coupling strength of each quantum term X. Every now and then in the 
sequence of consecutive displacement transformation the term X has the strongest coupling, and will be transformed 
away. By the arguments presented above, every next time we transform with X it will have a smaller coupling 
constant. 


CONVERGENCE OF THE METHOD 

As we have shown, the value of each coupling constant reduces through the consecutive application of displacement 
transformation. We will now estimate the speed at which the method converges, inspired by the works of D. Fisher[Tl[2]. 

For simplicity, let us focus on n = 4 order terms first, for a system with L sites. The Hamiltonian is, at each step 
of the renormalization procedure, given by 

ff(t) — ^ ^ 

a OLfB 


where the coupling constants are changing with every transformation. The total number of independent quantum 
terms are 


Nq = ^{L + 1)L{L-1){L-2). 


(24) 


Additionally, there are ^L{L — 1) classical terms at order n = 4, and L classical terms at order n = 2 that won’t 
change under renormalization. 

In the spirit of strong disorder renormalization group theory (SDRG), we define a probability distribution Pi{X) 
for the absolute value of the quantum interaction coupling constants X = \ \. Since we are explicitly looking at 

finite size systems, we must be careful about the interpretation of a probability distribution. A finite size system can 
be viewed as a specific realization with Nq terms taken from the distribution Pi{X). Given a continuous Pi{X), the 
expectation value for largest coupling constant becomes 

poo / pX 

T^ = Nq j dXX P,{X) f J dYP, 

For example, if Pi{X) is a uniform distribution between [0,1], the expectation value for T^ = 1 — + 0{L~^). 

Similarly, the expectation value for the second-to-largest coupling constant is 

dYP,{Y)\ dYP,{Y)^. (26) 


pOO 

n = NQiNQ-l) dXXP^iX) 

Jo 


No-1 


(Y) 


(25) 


The uniform distribution yields T' = 1 — which is the expected result. 

Each displacement transformation changes the distribution Pi{X). For a finite size system, we must pick the largest 
realized coupling constant and apply the transformation. This coupling constant becomes zero. The 

new distribution, by virtue of the fact that T^ was the largest realized coupling constant, should have zero weight for 
X >Ti. We obtain the distribution of the remaining couplings, 

^ P,{Y)dY^5{X). 


P,{X) = P,{X) - 0(X - T{)P,{X) + 


(27) 


The next step is to see how the given displacement transformation acts on the remaining terms. Only terms with 
maximal overlap will be able to generate new terms at the same order. Given any interaction there are Noveriap = 
(L — 4)(L — 5)/2 terms that have maximal overlap. As discussed in Sec. [j the overlapping terms always come in pairs, 
say Y and Z. The old interaction strengths Vy and are transformed into Vy cos X + Vy sin A and cos A — bV sin A. 
Consequently, a fraction Noveriap/Nq is changed by the transformation, and the new distribution function of couplings 
becomes 


P,+i(A) = P,(A) 1- 


X^overla 


N, 


Q 


+ 


N. 


overlap 


JO ^ ^ cos A,—/ ^sin A,- ^ 


N, 


Q 


cos A,: sin A,- 


(28) 


where Xi is the parameter associated with the displacement transformation. 

The remaining Nq — Noveriap — 1 terms are unchanged, hence Eqn. (27) and Eqn. (28) represent the change in the 
distribution function due to one displacement transformation. 

The method of consecutive displacement transformations converges if the expectation value of the largest coupling 
constant under the new distribution Pi^i{X) reduces significantly. We will now show that if initially only a sparse set 
of interactions have ’large’ coupling constants, the largest coupling constant reduces exponentially. If, on the other 
hand, the system has initially a large set of large coupling constants, the distribution PJX) will naturally become 
more and more sparse close to its upper bound. Once a sparse distribution is reached, the first argument again applies 
and a exponential decay sets in. 


Initially sparse set of interactions 


Consider a system where, initially, only 0{L‘^) of the quantum terms have a nonzero coupling constant. A typical 
initial distribution would be 


Po(A) = (i - - ^)- 

In the large L limit, the expectation value for the largest coupling constant is 

Tn = An 1 I — —IT 


L2 


(29) 


(30) 


Let us now apply Eqns. (27) and (28) to this distribution. The first step amounts to removing the largest coupling 
constant. We verify explicitly that we moved one coupling constant {0{L~'^) of the total number of terms) to zero 
strength, 


= 1 1 - L + i(x) + ^ (i - ^ I e(r„ - .Y). 


(31) 


The next step r equ ires knowledge of the transformation parameter Aq, which enters into the convolution of Pq with 
itself, see Eqn. (28). The self-convolution of Pq{X) gives for A > 0 explicitly 


1 


cos Ao sin Aq 


dZ Pn 


^X-Z\ ~ / Z 

ysin Ao, 


cos Aq / 

^o(l ~ ^o) 0(ro sin Aq — X) ^ ©(Tp cos Aq 

cos Ao 


To 


Tg sin Ao cos Aq 



0 < A < To sin Aq; 

To sin Ao < A < To cos Ao; 
cos Ao) — A, To cos Ao < A < ro(cos Ap -I- sin Ap); 
else 


(32) 


with Aq = — -p . The terms have a simple physical explanation. The first term with sin Ap is obtained because new 

terms are generated that had originally E = 0, but due to the transformation have become nonzero. The next term 
consists of terms that were originally present, but due to the transformation are reduced by a factor cosAp < l/v^- 
In both cases, it affects only a relative 0{L~‘^) number of terms. The final term on the last line consists of the pairs 
as discussed in Sec. [1 where both Vy and Vz were nonzero prior to the transformation. 
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This last part is dangerous, since it generates a weight of the distribution for X larger than Tg. However, even in 
the worst case scenario of a many-body resonance, which yields Ag = 7r/4, the integrated probability to find a new 
coupling constant larger than Tg is 


AIN, 


overlap / nz ^\2 0.68 


(33) 


We are saved by the sparseness of the initial set! In the new distribution Pi{X) the expectation value for the largest 
coupling constant is, in the large L limit. 


Ti^Tg 



(34) 


Since the top of the distribution remains 0{L ^) with our method, this estimate can be generalized to every step. 
Consequently, the coupling constant after i steps is given by 


r. ~ (^1 - ~ exp (-^/T') . 


(35) 


We have shown that for initially sparse interactions, the largest coupling constant will be reduced exponentially. This 
is consistent with, for example, the inset of Fig. |^of this Supplementary Information. 

Note that Eqn. (351 implies that the number of transformations required scales in a polynomial fashion with system 


size. The overall computational complexity of diagonalizing a Hamiltonian using displacement transformations, at a 
given order, is therefore also polynomial. 


Initially dense set of interactions 


The above arguments do not apply when our initial distribution has a dense set of large interactions. A typical 
example is the uniform distribution, 


Po(A) = ^0(Ao-A) 


(36) 


The expectation value for the largest coupling constant is now, ro = Ag(l—From the uniform distribution it is 
clear that there will be a macroscopic number of quantum terms that will, through the transformation, get a coupling 
constant larger than Fg. As Eqn. (321 shows, the new distribution function has a cut-off at Ai = rg(cosAg -l-sinAg) 
and it goes to zero in a linear fashion, Pi{X) ^ (Ai — X). 

The self-convolution of Pi will yield a distribution that will go to zero as P 2 {N) ~ (A 2 — A)^. In general, if Pi{X) 
goes to zero as (A^ — X)^, the next distribution has a tail of the form Pi^i{X) ^ (A^+i — The central limit 

theorem tells us that quickly the distribution will approach a normal distribution. 

A normal distribution is a prime example of a distribution with a sparse tail. Once the distribution Pi{X) has 
reached its Gaussian shape, we can revert to the arguments of the previous subsection to claim an exponential decay 
of the largest coupling constant. 

Notice that this argument is very similar to Fishers a posteriori justification for SDRG by introducing the infinite 
disorder fixed point [U [2]. Here the displacement transformation method seems to fail when there are many ’large’ 
interactions, however, the consecutive transformations make the tail of the coupling constant distribution so small 
that we reach a regime of exponentially sparse ’large’ coupling constants. 


Computational complexity 

We have just shown that for order n = 4 the largest coupling constant falls off exponentially with prefactor L^. 
For higher order interactions we observe that the relative density of terms with maximal overlap, Ngyeriap/NQ, scales 
as Consequently, to diagonalize a Hamiltonian up to order n given a fixed numerical precision, one needs 

to perform displacement transformations. Because there are L" quantum terms up to order n, the total 

computational complexity scales as 


CPU time~ 


(37) 
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FIG. 3: The computational time required to diagonalize the model of the main text up to order n = 4. We find polynomial 
dependence on the system size, with CPU time ~ 


We computed the actual CPU time needed to diagonalize the model introduced in the main text, up to order n = 4, 
as a function of different system sizes. The result is shown in Fig. which indeed satisfies a polynomial scaling in 
system size with the correct power L®. 


SIMPLE TEST-MODEL 


In the main manuscript we study the Anderson insulator with nearest neighbor interactions. Here we briefly discuss 
another, simpler model that we tested our method for. Consider a periodic chain of N sites with spinless fermions, 
with a random chemical potential on each site chosen uniformly between —W and W. The interactions couple four 
neighboring sites, with uniform strength V. The Hamiltonian is 

H = ^ ^ {c\ci+ic\^^Ci+3, + h.c.) . (38) 

i i 


One expects a localization-delocalization transition as a function of disorder strength W/V. 

We used a different code than used in the main manuscript. At each step we pick the quantum term with the 
largest coupling constant, and transforms it away. We neglect coupling constants smaller than numerical accuracy, 
set at e = 10“^^ in units where U = 1. To speed up the computation, each step we throw away all terms of order 
6 and higher. With this procedure, we indeed find that the magnitude of the largest coupling decreases rapidly, as 
shown in the inset of Fig. 

After order N'^ iterations, we have realized the classical Hamiltonian H = + VijUiUj. Within the model 


Eqn. (38), only next-nearest neighbor interactions are generated, Uj 5\i-j\=2- In this model, the structure of Uy is 


therefore not very enlightening to study the MBL phase. 

A better measure of the localization is to directly probe the locality of the new integrals of motion. To do so, we 
start out with a density operator Ui on a site, and transform it using the same transformations that diagonalized the 
Hamiltonian. The results will be of the form, up to quartic order. 


= C/In,C/ = ' 


jklm 


^jklm 


CnCkCi C-D 


(39) 


where U is the product of all the displacement transformations. Note that at quadratic order, rf = rij, consistent 
with the fact that the single-particle spectrum is unchanged by the interactions. 

Now there are various methods of determining whether a rf is quasi-local. Ref. [3] suggested to use the infinite- 
temperature overlap between lOMs at different sites. 


My = 4Tr(r/r/) - 1 (40) 

which is shown in Fig. for a = 12 chain for two different disorder strengths. One sees a clear indication of 
localization in the case of strong disorder W, and delocalization in the case of weak disorder. 
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FIG. 4: The normalized spread of the local integrals of motion for a = 36 chain with WjV = 5 and WjV = 0.5. The shaded 
area represents the standard deviation when averaging over all integrals of motion. The curve is normalized, so that the area 
under the curve equals one. Inset: The magnitude of the strongest coupling as a function of the number of transformations. 
It decreases exponentially, though much slower in the delocalized phase than in the localized phase. 



distance \i — j\ 


FIG. 5: The correlation function Mij — 4Tr(rfTj) — 1 as a function of distance, averaged over disorder realizations. This 
function expresses the localization of the local integrals of motion for a = 12 chain with WjV = 5 and WjV = 0.1. The 
shaded area represents the standard deviation. 


Another method, which is less time-consuming as it does not involve any trace, directly sums for each distance d 
the absolute value of the prefactors jal of terms that act on sites at distance d from each other. We computed this 
lOM spread on a A^ = 36 length chain, which is larger than state-of-the-art Exact Diagonalization studies can reach. 
The results of is shown in Fig. for two values of the disorder strength. Indeed, for strong disorder case we find 
localization whereas for weak disorder the lOMs have weight throughout the full length of the chain. 
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